Raster resampling simply refers to making changes to the pixel resolution of the raster. The term “resampling” implies that the pixel values are “sampled” and reassigned to the pixels at the new resolution. This operations often involves an interpolation method (nearest neighbor, bilinear, spline, min, max, mode, average etc). We will try three important functions for changing the resolution of a SpatRaster:
terra::resample - resample to match the resolution of another raster
terra::aggregate - resample from fine to coarse resolution
terra::disagg - resample from coarse to fine resolution
Raster resampling can be a critical step during multivariate analysis, where raster pixels must overlap to ensure the datasets from each raster corresponds to the same spatial domain. The following schematic helps illustrate the use of these functions:
We have so far worked with global surface soil moisture raster (from SMAP) and we are familiar with its spatial distribution across the globe. Now, let us consider this question:
How does the climate impact the spatial distribution of global soil moisture?
For this assessment, we will use Aridity index (AI), which is a popular approach for climate classification based on the relative availability of Mean Annual Precipitation (P) compared to the Mean Annual Reference Evapotranspiration (ET0) of a location. Aridity Index is defined as: \(AI= P / ET_0\). Here, ET0 is the maximum potential amount of moisture a hydrologic system can transfer to atmosphere through evaporation and transpiration. The value of \(P / ET_0\) progressively increases from arid to humid regions. Since P and ET0 can never be negative, AI is always >0. We have access to global raster of aridity index estimates (i.e. P/ET0) provided by (Zomer, Xu, and Trabucco 2022), at 5X5 KM spatial resolution.
One approach of ansering the question can be to compare pixel values of SMAP soil moisture with corresponding values of AI to establish the bivariate AI–SM relationship. So, we must first change the resolution of AI raster to match that of the SMAP soil moisture. For this, we will use raster resampling.
We will first import the raster and remove the negative fill values. We also plot a histogram of the raster to understand the probability distribution of the global AI values. Do we agree that most pixels in the raster range within 0-5?
library(terra)
AI=terra::rast("./SampleData-master/raster_files/P_over_ET0.tif")
terra::plot(AI)
# As we see AI raster has fill values. Prior to analysis, we remove the negative fill value
AI[AI<0]=NA
hist(AI, breaks=100, xlim=c(0,10)) 
# We observe most values in the raster are <4. So, we set the range of the plot accordingly
terra::plot(AI, main="AI=P/ET0", range=c(0,4))
First we will try an example for the terra::resample function to resample the pixels in the AI raster to match the spatial resolution of sm pixels using bilinear interpolation method. This method estimate new cell values as a weighted-average values of the adjoining pixels. The weights are calculated using the distance of the centriod of the target pixel from the adjoining cells. In addition to the bilinear approach, terra::resample has several interpolation methods available as options such as near (for nearest neighbor interpolation), cubic, sum, min, max, average , rms (root-mean square) etc.
Once AI raster is resamapled to match sm, we would then be able to generate scatter plots between the two rasters, and evaluate the relationship between aridity index (climate) and soil moisture. We see that the soil moisture increases as AI increases before reaching an asymptotic value as AI>0.65 (humid climate, indicated by red vertical line in the plot). This relationship follows the famous Budyko formulation of energy and water limits on terrestrial water balance (Chen and Sivapalan 2020). For illustration, we use a simple formulation of Budyko curve to represent the non-linear interrelationship between AI and soil moisture given as:
\[ SM=AI/(1+AI)^{1.5} \]
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
aridityResamp=terra::resample(AI, # Raster to be resmapled
sm, # Target resolution raster
method='bilinear') # bilinear interpolation method
plot(aridityResamp,sm ,
xlim=c(0,3), ylim=c(0,0.6),
xlab="P/ET0", ylab="Soil Moisture",
pch=19)
curve(x/((1+x)^1.5), col="blue",lwd=3, add = T)
abline(v=0.65, col="red", lwd=3) 
Food for thought: Can we do this analysis had we used
AIraster instead ofaridityResamp? What will be the output is we replacearidityResampwithAIin the code above.
In the AI plot above, we see the expected climate pattern for terrestrial landmass. The equatorial regions receive abundant precipitation, for example the Amazon forest and Southeast Asia, and have high AI values. Regions like Sahara desert, Central Australia, Western US are have hot and dry climate, and show low values of AI.
Now that we have seen the application of terra::resample function, let us now try terra::aggregate and terra::disagg when we know the factor of (dis)aggregation to be used in each direction. Several functions are available for raster aggregation including mean, max, min, median, sum, modal, sd. Disaggregation can use either near or bilinear methods.
library(terra)
# Import SMAP soil moisture raster
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
# Original resolution of raster for reference
res(sm)[1] 0.373444 0.373444#~~ Aggregate raster to coarser resolution
SMcoarse = terra::aggregate(sm, # Soil moisture raster
fact = 10, # Aggregate by x 10
fun = mean) # Function used to aggregate values. We use within-pixel mean
res(SMcoarse)[1] 3.73444 3.73444plot(SMcoarse, main="Aggregated soil moisture")
#~~ Disaggregate raster to finer resolution
SMfine = terra::disagg(sm,
fact=3,
method='bilinear')
res(SMfine)[1] 0.1244813 0.1244813plot(SMfine, main="Disggregated soil moisture")
United Nations Environment Program (UNEP, (Nash 1999)) divides global climate in five classes based on AI as below:
Several applications require discrete classification of variables. A discrete attribute table can be added to a raster which serves as a reference point for discrete classification of the continuous variable in the raster. As an example, we will follow class breaks and names given in Table 1 to add climate attribute to the AI raster.
# Import AI raster and remove negative fill value
AI=terra::rast("./SampleData-master/raster_files/P_over_ET0.tif")
AI[AI<0]=NA
# Breaks for each climate class (Table 1)
class_brk= c(0, 0.03, 0.2, 0.5, 0.65, 10)
# Labels for each climate class
class_names=c("Hyper arid", "Arid", "Semi arid",
"Sub humid", "Humid")
# Divide the cell values in the AI raster into distinct levels
attributes=cut(values(AI), # Notice that we apply cut after extracting raster values.
breaks = class_brk,
labels =class_names )
# Add attributes to the SpatRaster as climate class
AI$climate = attributes
plot(AI$climate, plg = list(loc = "bottomleft")) 
Let’s explore using a spatial polygon/shapefile for summarizing a raster (in this case, global SMAP soil moisture) by using extract function from the terra library. We will also transform global aridity raster to a polygon using the function as.polygons to find the mean soil moisture values for each aridity class.
First, we will use the IPCC shapefile to summarize the soil moisture raster.
library(terra)
# Import SMAP soil moisture raster from the downloaded folder
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
# Import the shapefile of global IPCC climate reference regions (only for land)
IPCC_shp = vect("./SampleData-master/CMIP_land/CMIP_land.shp")
#~~~ Using shapefile to summarize a raster
sm_IPCC_df=terra::extract(sm, # Spatraster to be summarized
IPCC_shp, # Shapefile/ polygon to summarize the raster
fun=mean, # Desired statistic: mean, sum, min and max
na.rm=TRUE) # Ignore NA values? TRUE=yes!
head(sm_IPCC_df) ID SMAP_SM
1 1 0.2545766
2 2 0.3839087
3 3 0.2345457
4 4 0.3788003
5 5 0.2383580
6 6 0.1475145#~~~ Extract cell values for each region
sm_IPCC_list=terra::extract(sm, # Raster to be summarized
IPCC_shp, # Shapefile/ polygon to summarize the raster
fun=NULL, # fun=NULL will output cell values within each region
na.rm=TRUE) # Ignore NA values? yes!
# Apply function on cell values for each region
sm_IPCC_mean=lapply(sm_IPCC_list,mean) # Returns a list of regional means
sm_IPCC_mean=purrr::map(sm_IPCC_list,~ mean(.x, na.rm=TRUE)) # Returns a list of regional means
#~~ Try user defined function
myfun=function (y){return(mean(y, na.rm=TRUE))} # User defined function for calculating means
#~ Implement function using lapply and map
library(purrr)
sm_IPCC_mean=lapply(sm_IPCC_list,myfun) # Returns a list of regional means
sm_IPCC_mean=purrr::map(sm_IPCC_list,~ myfun(.x)) # Returns a list of regional means
sm_IPCC_mean=unlist(sm_IPCC_mean) # Unlist to return a vector
head(sm_IPCC_mean) # Is this the same as the previous result? ID SMAP_SM
21.1513647 0.2088026 In the next example, we will convert global aridity raster into a polygon based on aridity classification using terra::as.polygons function.
Global aridity raster has 5 classes with 5 indicating humid and 1 indicating hyper-arid climate. We will use this polygon to extract values from the SpatRaster and summarize soil moisture for each aridity class.
#~~~ Convert a raster to a shapefile
aridity=rast("./SampleData-master/raster_files/aridity_36km.tif") #Global aridity
# Convert raster to shapefile
arid_poly=as.polygons(AI$climate) # Convert SpatRaster to polygon and then to sf
# View aridity polygon
library(mapview)
mapview(arid_poly)